REPORT  No.  907 


c«c0( 


On  The  Choice  Of  Mesfir-”S-d; 
In  The  Integration 
Of  Ordinary  Differential  Equations 

Approved  for  public  release;  distribution  unlimited. 


BORIS  GARFINKEL 


TECHNICAL  LIBRARY 
DRXBR-LB  (Bldg.  305) 

ABERDEEN  PRO/ING  GROUND.  HD.  21005 


DEPARTMENT  OF  THE  ARMY  PROJECT  No;  503-06-002 
ORDNANCE  RESEARCH  AND  DEVELOPMENT  PROJECT  No.  TB3-OQ07K 

BALLISTIC  RESEARCH  LABORATORIES 


r*' 


ABERDEEN  PROVING  GROUND.  MARYLAND 


bo  "Initial  distribution  has  been  made  of  this  report  in  accordance 
with  the  distribution  list  contained  herein o Additional  distribution  with- 
out recourse  to  the  Ordnance  Office  may  be  made  to  United  State's  military 
organizations s and  to  such  of  their  contractors  as  they  certify  to  be 
cleared  to  receive  this  report  and  to  need  it  in  the  furtherance  of  a 
military  contract o" 


DEPARTMENT  OF  THE  ARMY 
U.S.  ARMY  ARMAMENT  RESEARCH  AND  DEVELOPMENT  COMMAND 
U.  S.  ARMY  BALLISTIC  RESEARCH  LABORATORY 
ABERDEEN  PROVING  GROUND,  MARYLAND  2100S 

REPLY  TO 
ATTENTION  OP 

DRSMC-BLS [A)  22  September  1983 

SUBJECT:  BRL  Report  No.  907  s 


SEE  DISTRIBUTION 


1.  Reference  BRL  Report  No.  907,  AD  number  044458,  "On  the  Choice  of  Mesh 
in  the  Integration  of  Ordinary  Differential  Equations",  by  B.  Garfinkel,  dtd 
April  1954,  UNCLASSIFIED. 

2.  Request  that  you  mark  your  copy  of  the  subject  report  with  the  following 
distribution  statement; 

Approved' for  public  release,  distribution  unlimited 


and  attach  a copy  of  this  letter  to  the  report. 


Chief 

Security  Office 


BALLISTIC  RESEARCH  LABORATORIES 

REPORT  HO.  907 
April  195U 


ON  THE  CHOICE  OF  MESH  IN  THE  INTEGRATION  * 
OF  ORDINARY  DIFFERENTIAL  EQUATIONS 


Bori a Garflnkel 


TECHNICAL  LIBRARY 
DRXBR-LB  {Bldg.  305) 

ABERDEEN  PROVING  GROUND,  HD.  21005 


Department  of  the  Amy  Project  No.  503-06-002 
Ordnance  Research  and  Development  Project  No.  TB3-0007K 


ABERDEEN  PROVING  GROUND,  MARYLAND 


BALLISTIC  RESEARCH  LABORATORIES 


REPORT  HO*  907 


BG  arfinkel/ddh 

Aberdeen  Proving  Ground,  Md„ 

April  1 95k 


x ON  THE  CHOICE  OF  MESH  IN  THE  INTEGRATION 
OF  ORDINARY  DIFFERENTIAL  EQUATIONS 


ABSTRACT 

The  optimum  mesh  at  point  x is  defined  here  as  the  interval,  h(x), 
which  minimizes  the  time  of  integration  of  an  ordinary  differential 
equation  while  maintaining  a prescribed  bound  of  the  error®  An  attack 
upon  the  problem  of  constructing  such  a mesh  is  carried  out  here  with 
the  aid  of  the  Calculus  of  Variations*  The  Euler  equations  are  derived 
and  discussed,  and  the  results  are  extended  to  a system  of  differential 
equations* 
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Intentionally  left  blank. 


IHTROEU  CTION 


In  the  numerical  integration  of  ordinary  differential  equations  by 
a sfcep-by-7fcep  method  it  is  sometimes  advantageous  to  change  the  steps 
h,  as  the  integration  proceeds.  This  suggests  the  following  problems 
“Given  the  equation 

y1  “ f(x»y)»  y(o)  y0,  o (1) 

and  a numerical  method  whose  local  truncation  error  is  of  order  hk,  k - 2, 
required  the  function  h(x),  such  that  1)  the  accumulated  error,  6y, 
does  not  exceed  a prescribed  bound,  E-fc,  and  2)  the  time  of  integration 
is  minimzed% 

Provide^  6y  is  sufficiently  small,  it  will  satisfy  the  variational 
equation 


4 Sy  = fy6y  + V/h,  (2) 

whore  is  the  local  error  per  step.  If  £ is  an  upper  bound  of  the 
rate  of  the  local  error;  iBe«, 

|r)/4^J  (3) 

the?}  5 y has  an  upper  bound,  E(x),  satisfying  the  differential  equation1^ 


3*  - f^E  -6-  0,  E(0)  = 0, 


(h) 


whose  solution  is 


E(x)  = 


x 

n 


x 

exp  J f (a,y(a))  ds€(£ ) d£.($) 


bound  £ - of  the  error  rate  can  be  estimated  as 
+ j 

€ z - R/h  + hk_1^k  + hk^k+1  + .,,, 


(6) 


VJhere  R is  a bound  of  the  local  rounding  error  and y,  (x)sp,  _(x)  are 

” K , K+J_ 

functions  of  the  zeroth  order  in  h*  If 


R ^ (k-l)ykk+1/  V k+1 


it  is  possible  to. choose  h such  that 


h 5 
m - 


R/(k-i)  yk|1A<h«  V^k +r 


(7) 

(3) 
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Then  the  second  term  in  (6)  dominates  the  third , and  € khk”^  so 

that  a new  bound  can  be  conveniently  chosen  as  £ «*  kh'  , the  subscript 
k of  V having  been  dropped*  If  the  rounding  error  is  negligible,  a more 
practical  bound  would  be 

€ = hk"1  y , (9) 

which  we  shall  adopt  here*  Dn  the  other  hand,,  if  the  rounding  error  is 
not  negligible,  the  results  to  be  derived  wili  remain  valid  provided  Y 
is  replaced  by  k \jJ  wherever  it  occurs* 

The  simple  problem  of  minimizing  E(x),  regardless  of  the  time 
consumed,  is  seen  from  (£)  to  be  equivalent  to  that  of  minimizing 

“ R/h  + h^”ly'0  Two  properties  of  the  solution,  furnished  by  the 

function  h =*  h (x),  may  be  noted  in  passings  1)  €-(h  ) 0 khk”^y,  2)  Rs 
in  1 m m 

hk\y  = (k-l)s  lj  i,e*  the  local  error  bound  is  partitioned  between  the 

m 2 ) 

rounding  and  the  truncation  in  the  ratio  (k  - 1)}  1,  The  latter  result  - 

should  be  contrasted  with  the  popular  notion  that  this  ratio  should  not 

exceed  unity® 

The  function V (x ) is  listed  below  for  some  of  the  more  commonly 
used  methods  of  step-by-step  integration* 


Table  1 


Method 

k 

X 

Remarks 

Euler 

2 

\ y,fl 

' ( 

Modified  Euler ^ 

3 

il  r"  -sV'" 

one  iteration 

tf  ft 

3 

two  iterations 

Kutta 

h 

Runge-Kutta^ 

1 ...  

5 

c(x)y^^ 

The  modified  Euler  method  with  one  iteration,  called  by  some  the  Heun 
method,  can  be  cited  as  an  exception  to  the  general  rule  that  y;  is  of  the 
farm 

V (x)  * c(x)y^  (x),  (10) 

where  y^k^  is  the  derivative  of  order  k and  c(x)  depends  on  f(x,y)  and 
the  method  used® 
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If  new  variables  H(x),  S(x)  are  now  defined  as 


X 

H(x)  - exp  (-  0 

f (^)d Z)y)(x)Mo) 

tf 

(11) 

S(x)  3 exp  (- 

fv(£)  d£  ) S(x), 

1 

(12) 

(lj)  can  be  rewritten 

0 

S *-^(0)hk"1  H3*  = 

0,  S(0)  = 0. 

(13) 

The  time  t * of  integration  for  a given  method  is  proportional  to  the 
number  of  integration  steps'.  Therefore, 

x 

t(x)  = J d£/h(£),  (l/i) 


or 


(ISO 


the  three  unknown  functions  h,  S,  and  t are  thus  connected  by  two  non- 
bolonomic  constraints  (13)  and  (l5)«  The  system  therefore  has  one 
degree  of  freedom,  which  can  be  realized  by  an  arbitrary  choice  of  h(x). 


VARIATIONAL  APPROACH 

We  identify  our  problem  with  the  Problem  of  Belza^  in  the  Calculus 
of  Variations:  "Required  the  function  h(x)  in  the  prescribed  range 
(x^  = 0,  Xg  = X ) s which  minimizes  the  function  g t(xg)  and  satisfies 

the  differential  constraints 


- V (0)hk-1Hk(x)  - 5«  = 0, 

* h”1  - t«  =0, 

subject  to  the  end-conditions 

S(0)  «t(0)  = 0,  S(X)  = S*  =E*expC“ 


X 


(16) 


(17) 


The  Lagrangian, 


must  satisfy  the  Euler  equations 


3x  Fz'  * Fz.  5 i “ 2,  3j 

a i 

(19) 

srith 

z^  ~ h,  z2  = S,  *r  t. 

(20) 

Then  (l?) 

becomes 

\1  (k-1 ) y ( 0 -\2h"2  * 0, 

X^  = const.. 

(21) 

^2  * const.. 


leading  immediately  to 
h = h(0)H 


-1 


(22) 


h(0) 


X2/(k-l)  1|J  (0)^ 


— = const , 

k 


The  parameter  h(0)  is  determined  from  the  end-conditions  as 


(23) 


and  t(X)  from  (22.1),  (l6.2)  as 

X 

t(x) » r hco  /h(o). 

J0  - .:‘- 

Xt  can  be  shown  that  t z T,  where  T corresponds  to  a constant 
interval,7i.  For  we  have  from  (l6) 

x , 1 


(2U) 


S*Y  V (o)  J ^ (£)d£ 


l7=I 


(25) 


t = xa; 
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6) 

Hence,  with  the  and  of  the  H&lder  inequality  ' for  k-^l,  there  follows 

(26) 

whore  H and  are  the  mean  values  of  H and  respectively,  on  the 
interval  (0,X)„  This  result,  of  course,  would  be  expected  if  our 
solution  furnishes  the  absolute  minimum  of  t. 

From  (l6)  and  (22,1)  can  be  deduced 

dS/dt  « consti  (27) 

x 

E(x)  =■  t exp  p f y (£)d£const,, 

J0 

which  admits  of  the  following  interpretation;  "AH  points  of  an  optimum 
mesh  make  equal  contributions  to  the  final  truncation  error", 

THE  EFFECT  OF  DISCONTINUITIES 

Since  the  use  of  a continuous  h(x)  is  physically  impossible,  it  may 


be  replaced  by  some  convenient  step-function,  such  as 

\(x)  = (0)2m,  m = 0,  + 1,  + 2,...,  (28) 

with  the  additional  constraint 

5/3  = m'  = 0,  (29) 

the  augmented  Lagrangian, 

F = ^ (y  (0)hk_1Hk-S« ) + \2(h-1-  t*  ) + X^m*  (30) 

7) 

must  satisfy  the  Comer  Condition  , 

AFZ,  = 0,  A(F  - zJ  Fz,)  - 0,  (31) 

i i 

where  A denotes  a "jump"j  e0g,, 

AF  = F+  •»  F_.  (32) 

From  (31)  and  (30)  is  obtained 

A =AX2  = AX 3 = 0 

A (Xx  VCO)^”1  + X.h'1)  *=  0,  (33) 
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If  f is  of  class  C 


k-1 


then 


A f - 0,  A H = 0, 


(3U) 


and  we  deduce,  with  the  aid  of  (28), 

A / A2‘m  » const,  H"lc  (3$) 

Finally,  if  m(0)  - 0 and|^m|=  1,  (3£)  yields  I^2‘m  = 1 at  comers,  and 

r. 


ra(x)  = - log  H(x)/log  2 


(36) 


holding  everywhere,  the  symbol  [ denoting  the  integral  part  of  a number 
h(0)  and  t(l)  can  now  be  calculated  from  (17)  as  1 

x Ira 


Mo) 


7y(o)  J 


,(k-l)m(  £,)  „k(  £)d  6 


H 


(37) 


>-*»(£) 


d £ A*  (0). 


THE  EFFECT  OF  INEQUAUTTES 

Since  R,  ^k+l  occ’urrinS  (9)  are  generally  difficult  to 

calculate,  some  other  convenient  bounds  may  be  imposed  on  h in  the  general 
form 

a^h  (x)  —bo  (33) 

Furthermore,  instead  of  prescribing  the  terminal  error  bound,  E(X),  there 
may  be  imposed  a bound  on  the  function  E(x): 

E(x)4E*.  (39) 

Inasmuch  as  h(x)  is  discontinuous,  the  usual  Tangency  Condition  does  not 
hold,  so  that  the  inequalities  (38),  (39)  must  be  enforced  directly. 

It  is  to  be  noted  that  putting  an  upper  bound  on  h prevents  the  break- 
down of  the  solution  when  if*  (x)  = 0,  leading  to  H(k)  * 0 in  (H)  and  h(x)  =°° 
in  (22). 

THE  SUFFICIENCY  CONDCTION 

For  a strong  relative  minimum  the  Sufficiency  Conditions^  are  I,  II* 

HI*  , IV,  described  below*  The  Multiplier  Rule,  I,  is  satisfied  by  the 
solution  of  the  Euler  equation  provided  and  7^  are  not  both  zero  and 

the  Transversality  Condition  is  satisfied".  The  latter,  in  our  problem 
reduces  to 
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«t  * Ff 


3 


m 


= 0 


or 


X2(x2)  - 1.  CUl) 

The  Weierstrauss  Condition,  II* , expressed  in  terms  of  the  E- 
function,  is 

£(z8 ,2* ) - FCx,z,Z«)  - F (x,z,z«)  - (2»-z»  )Fz,(xsz,z*  )^0,  (1|2) 

where  the  slope  functions  z*  are  functions  of  the  field  coordinates 
SjZ  in  the  neighborhood  of  the  solution  and  2*  ^ zf  * If  the  derivatives 
of  certain  variables  z,  are  lacking  in  F,  then  such  z.  must  be  classed 
a slope  function  and  K(h2)  modified  accordingly,,  Inour  problem,  there- 
fore, the  slope  functions  are  h,  5* , t',  and  (1*2)  becomes,  in  view  of 
(30),  (22) 

E » X-^  (OjH*  + Xjdh-1 


“ Tg-TJKtln^T'  L(1  + ” 1 “ k P 0*3) 

where  6h  denotes  a strong  variation  and 

P =6h/h.  (Wi) 

we  observe  that  h-^0,  p^-l,  k^2,  and 

(1  + p)k  - 1 - kp^O  if  - l*p  J 0,  k^O  (1*5) 

Thus  (1*2)  is  equivalent  to  the  requirement 

X2(x)  0,  (1*6) 

which  is  obviously  satisfied,  since 

X2(x)  = 1 (1*7) 

in  virtue  of  (1*1)  and  (21) 6 

The  Clebsh  Condition,  IH 1 s is 

i6z*^°  «*»> 


for  all  5zJ^  satisfying  the  differentiated  equations  of  constraint.  In 
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our  problem  h(x)  is  the  only  slope  function  contributing  to  the  quadratic 
form,  so  that  (U8 ) becomes 

Xgkh"5  6h2^*  0,  (U9) 

which  is  automatically  satisfied,  since 


1,  k ^2,  h^  0, 


The  positive-definiteness  of  the  second  variation,  which  is  required 
by  Condition  IV* , in  our  problem  reduces  to 


d‘ J = j \2kh~3  6h2dx  ^ 0, 
0 


(?o) 


since  x^  and  x2  are  fixed,  while  g and  the  end-functions  are  linear  in 

their  arguments*  Since  all  the  conditions  are  satisfied,  the  unique 
solution  (22)  furnishes  an  absolute  minimum  of  t"« 


AN  APPROXIMATE  SOLUTION 


,0c) 


If  f , yv  > c(x)  are  r placed  by  some  constant  bounds  defined  hy 

v 

means  of 


|/(x)  | £ ML1'”1,  |c(x)|^c, 

then  (10)  becomes 

V (x)  = cML15"1  ■ y (0) 

In  terms  of  the  dimensionless  quantities 

x<  * Lx,  h*  - Lh, 

(11),  (12)  reduce  to 

H - a“X/\  S - E(x)e“x 

cm  dropping  the  primes,  while  (22-25)  become,  respectively 


(5D 

(52) 

(53) 
(5U) 


h = h(0)e 


x/k 


h(0) 


LEfc 

Me 


-X 


1 

FT" 


k(l-e'X'^) 
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too 

K - 


k 

TO) 


(1  - e-*A 


(55) 


t(x)  = x/K  , 

for  the  discontinuous  solution  (28),  (36),  (37)  become,  respectively. 


h^(x) 


(0^(>:)  ^ 


m(x)  « £x/k  log  2^, 


h*(°) 


LE*e"X 


1 


Me  2(1-2"*)  (l-e"X^) 


(56) 


x _ k log  it  .-xA\ 

**  ~ ^TOT'  (1  " 6 }* 


A comparison  of  the  discontinuous  and  the  continuous  solutions  yields 

h»  (x):  h(x)  » 2"e  “x^  f(k),  (57) 

J t - f"1  (k)  log  b, 

where  f(k),  defined  by 


1 


f(k)  = k/2(l-2“k) 


FT 


(58) 


is  a weak  function  of  k*  It  has  the  properties 

f(0)  “ log  it,  f(l)  * e/2,  f(<»)  = 1,  f'^0  (59) 

and  is  tabulated  below. 


Table  2 


k 

r> 

1 

2 

3 

it 

5 

6 

f(k) 

1.39 

1.36 

1.33 

1.31 

1.29 

1.27 

1.25 

t*:t 

1.00 

1.02 

1.05 

1.06 

1.08 

1.09 

i.n 

In  the  range  2 £ k ^ 00 

1 6:  f (k)  it/3,  (60) 
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m •oc/ic 

while  2 e ' is  a periodic  piecewise  raonotonic  function  oscillating 
between  the  values  l/2  and  1 with  the  period  x » k log  2*  Consequently, 
h^:  h oscillates  in  the  range 

1/2  £ f(k)/2  4z  h#:  h*f(k)*U/3,  (61) 


so  that  the  continuous  and  the  discontinuous  solutions  are  interlaced,, 


If  X is  approximated  by  an  integral  multiple  of  the  period  k log  2, 
we  may  put 


X = N k log  2 
and  deduce  from  (55 ),  (56) 


t„x  1 


2(1-2"N) 
Jt  ' 


2(l-2“k)  (1-2^) 


1-2 


-Hk 


1 


■=  0(tf,k). 


(62) 


(63) 


to  be  noted: 


special  values  and  an  a 

0(1,10  “ 1, 

0(°°,k) 

0(N,k)  - | 1 

2(l-2'k) 

1 

ft 


if 


(6k) 


Consequently,  as  N ranges  from  1 to  » the  “relative  gain" , 1 " 0,  of 
the  optimum  mesh,  h (x ) v in  comparison  with  the  constant  mesh,  "H, 

ranges  from  0 to  1*  It  should  be  noted,  however,  that  the  use  of 
constant  bounds  in  (5l)  may  lead  to  a solution  h(x)  which  is  fantastically 
smaller  than  the  one  based  on  the  local  bound  £(x),  varying  from  point  to 
point o For  this  reason  the  results  of  this  section  must  be  used  with 
caution* 


ON  THE  CHOICE  07  METHOD  OF  INTEGRATION 


If  we  do  not  restrict  ourselves  to  a particular  method,  the  preceding 
results  enable  us  to  choose  the  optimum  k*  By  making  the  reasonable 
assumption  that  the  time  of  integration,  per  step  is  proportional  to  the 
order  of  accuracy,  k-1,  we  can  write  t as  an  explicit  function  of  k,  in  the 
form 


t = (k-1) 


Hu (S)d£  ft.  (o)/s 


-it- 


1 

ft 


(65) 


where  is  a known  function  of  x,  and  proceed  to  minimize  t with  respect 

Ik 


to  k.  For  example,  if  Hk(x),  yk(0)  should  be  constants  independent  of 
k,  the  optlimm  k = k^  is 

k*  »[l  + log  aJ  3 

1=  ^k(0)RkX^  (66) 

As  one  would  naturally  expect,  increases  with  the  range  X and  with 
the  precision  required,, 

In  iterative  procedures  the  number,  y (x),  of  cycles  per  stop  can 
also  be  optimised,,  This  can  be  achieved  by  introducing  v as  a factor 
of  t and  regarding  € (h,v)  as  a known  function  of  h,  v,  as  well  as  x® 

Then  h(x)  and  v(x)  are  determined  from  the  Euler  equations® 

heh  * v£v  “ 0 5 
X 

log  (h2  £ h/v)  f ( £ ) d$  + const®  (6?) 

0 * 

A SYSTEM  OF  DIFFSIEOTIAL  EQUATIONS 
The  construction,  of  optimum  mesh  can  be  extended  to  a system 


y±U)  - ^(x^)  s 7±(0)  =>  7ioi  (68) 

i , J = I3  » ®nB 

The  Lagrangian  can  now  be  written 

F * £(£+  JE  - E* ) + X (h-1-t»),  (6 9) 

where  J is  the  Jacobian  matrix 

Jij  * (70) 

E and  £ are  column  matrices  satisfying  the  relations 

6Fi-Si» 

(71) 

and  fi(x)  is  a row  of  Lagrange  multipliers*  Generally*  in  this  section  a 
bar  placed  above  a letter  will  denote  matrix  transposition.  Let  the 
end-conditions  be 

E(o)  a 0,  t(0)  * 0,  £(X)  I - E*  «*  0 (72) 


15 


The  Euler  equations  yield 


H1  “ - H 


h = h (0)  H"1 

H£|^(0)  Y(°)  | 1/k» 

The  solution  of  (73el)  is 

»U)  a p!(0)  2 (x)5 


(73) 


(7U) 


where  2 is  the  matrix  solution  of 


z>  «-n,  2(o)  = i. 


(75) 


X being  the  unit  matrix#  Since  J 
soon  as  y(x)  is  known*  it  remains 
constants  ^(0)*  h(0)*  satisfying 


and  Z became  known  functions  of  x as 
to- determine  the  n independent 
the  constraint 


(76) 


The  latter  condition  may  be  imposed  without  any  loss  of  generality  a since 
(7k)  is  a linear  homogeneous  equations  |i(0)  can  be  found  with  the  aid  of 
the  Tran sver  sail ty  condition  Y 


+ F. 


0S 


(77) 


where 


g + a 


i 


= t + a ( | E | -E% 


(78) 


a being  a constant  Lagrange  multiplier#  From  (77)>(78)*  (69)  follows 

n(x)  - a E(x)/  | e(x)  | , (79) 


which*  in  view  (7$)s  (77)  yields 

^l(o)  = tE(x)z“1(x)y/  EOO*-1^).  V (80) 

Although  E(X)  appearing  in  (80)  is  not  immediately  available 5 it  is 
possible  to  calculate  an  approximation  by  extrapolating  y(X)  to  "zero- 
mesh"#  An  iterative  process  may  then  be  set  up  between  (80)  and  (73)s 
which  may  converge  to  yield  the  value  of  p,(0)„  Finally * h(0)  may  be 
determined  from  the  end-conditions  j with  the  aid  of 

x 

E(x)  = Z^Cx)  J Z(£)V  ($,)  hk“1(^)  d£  (81) 

0 
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9 


and  of  (73  )s  as 


h(0)  = f Z(X)E(X)/( 

7 J,  |jiCO}Z(S)  f (£)  I1-1* 


I FT 


Obviously,  if  n = 1 


J - f , Z * exp  (. 


x 


1 


d£  ),  n(0)  - 1 


(82) 


(93) 


so  that  (81),  (82),  and  (73®3)  reduce  to  (5)  (23)  and  (ll),  respectively® 

NUMERICAL  EXAMPLE 

Required  the  numerical  solution  of  yB  - 2y  -e2",  v(0)  “ 1 by  the  method 
of  Runge-Kutta,  with  a precision  of  E#  ® 1 at  X =®  10.  tu  Since  the  formal 
solution  is  known  to  be  y « qk9  it  is  possible  to  write  down  the  theoretical 
optimum  mesh®  Here  (x)  = c(x)ex  , ^ (0)  = c(o)  From  Table  1, 

f ® 2 so  that  from  (ll)s  (12),  (36)3  (37)  there  follows? 


H = e"X/7c  3 


VO 


3*  - e“2x,  m - f x/k  log  2~j , 

■ 1 
S~ 


e(0)  2 (1-2”k ) (1-2“N) 


FT" 


and  from  (25) 


**  -k^t  (2k  ^ 2>  f1  - 2'h  )> 


Tt 


c(o)  d-2"r) 

X = x/5  * 


1 

FT 


(6U) 


(85) 


In  our  example  k = 5*  hence  M ® 3 ffom  (62)|  c(x)  = 0,073  ~ c(0)  from 
formula  (13)  of  reference  3)®  The  numerical  values  now  becomes 


m(x)  = 
h*(0)  = 


[ x/3,li7 
0*0093  , 


t = 65l  steps, 
"E  - 0,00.06  s 

1 = 981  steps. 


(86) 


giving  the  optimum  mesh  a relative  gain  of  3h%* 
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Of  course*  in  numerical  work  the  solution  y(x)  is  not  known  a 
priori®  The  optimum  mesh  can  be  constructed  by  a computing  machine,, 
using  a rough  integration  of  the  differential  equation  y*  = f {x*y) 
with  an  arbitrary  initial  h(0)*  combined  with  the  following  algorithm; 


x 


- 

- 

t 

Am  «■ 

( 

fy(  £ )d<£-  log  r\/r)  (o) 

) / k log  2 

0 


V (Ay(1)  -^y(2)/(l-2rk+1)s 


m(xi+i)  = mU^)  +Am* 
h = h(0) 


(87) 


*±♦1  * *i  + 


1 


Here  V)  is  an  estimate  of  the  local  truncation  error 3 which  is  calculated 
by  means  of  extrapolation  to  "zero  mesh"?).  AyW  andAy(^'are  the 
increments  of  y due  to  step  h and  a succession  of  two  half -steps,,  h/2a  *. 
. (87®1)  is  a consequence  of  (36)  and  (11).  The  constant  h(0)  satisfying 
the  condition  E(x)  = E*  can  be  found  by  some  suitable  Error  Control 
procedure^)® 


CONCLUSIONS 

We  have  shown  that  an  optimum  mesh  exists  and  have  described  the 
process  whereby  it  could  be  constructed®  Generally,  the  saving  of  time 
depends  on  the  differential  equation  to  be  integrated,  the  method  u3ed* 
and  the  range  of  integration®  As  would  be  expected*  the  saving  increases 
with  the  range*  If  the  latter  is  sufficiently  great*  the  gain  will 
exceed  the  cost  of  additional  labor  necessary  for  the  construction  of  the 
optimum  mesh*  and  will  fully  justify  its  use. 


BORIS  GAR5TNKEL 
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